% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

load('no_policy')

%% Generate random shocks

options = optimset;
T       = 100000;               % simulation periods
t_drop  = 1000;

rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);
dv   = zeros(T+1,1);
rkhats_next  = zeros(T+1,1);
rkhatss_next = zeros(T+1,1);
pd_next  = zeros(T+1,1);
pr_next  = zeros(T+1,1);

rkhat_insol  = zeros(T+1,1);

%% Simulation

zthre     = 0;
num_insol = 0;

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;
lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        % insolvency check
        rkhat_insol(t) = log ( rbv(t-1) * (1-1/lv(t-1)) );
        if rkhat(t) < rkhat_insol(t)
            num_insol = num_insol + 1;
        end
        
        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next(t)  = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next(t) - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;

        rkhatss_next(t) = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda           ) - (1-delta) );
        z_temp = (rkhatss_next(t) - erk(t)) / sigma_rk;
        pd_next(t) = normcdf(z_temp);
        pr_next(t) = prv(t) / pd_next(t);

        dv(t) = kv(t) - nv(t);

    end
    % End of T period simulation

avg_c_nopol  = mean(cv(t_drop+1:T))
avg_y_nopol  = mean(yv(t_drop+1:T))
avg_h_nopol  = mean(hv(t_drop+1:T))
avg_l_nopol  = mean(lv(t_drop+1:T))
avg_k_nopol  = mean(kv(t_drop+1:T))
avg_i_nopol  = mean(iv(t_drop+1:T))
avg_n_nopol  = mean(nv(t_drop+1:T))
avg_rb_nopol = mean(rbv(t_drop+1:T))
avg_pr_nopol = mean(prv(t_drop+1:T))
avg_pib_nopol= mean(pibv(t_drop+1:T))
avg_d_nopol  = mean(dv(t_drop+1:T))
avg_rkhat_nopol  = mean(rkhat(t_drop+1:T))
%avg_spread_nopol = mean(exp(rkhat(t_drop+1:T))-rbv(t_drop+1:T))
avg_spread_nopol = mean(exp(erk(t_drop+1:T))-rbv(t_drop+1:T))
avg_erk_nopol= mean(erk(t_drop+1:T))

std_c_nopol  = std(cv(t_drop+1:T))
std_y_nopol  = std(yv(t_drop+1:T))
std_h_nopol  = std(hv(t_drop+1:T))
std_l_nopol  = std(lv(t_drop+1:T))
std_k_nopol  = std(kv(t_drop+1:T))
std_i_nopol  = std(iv(t_drop+1:T))
std_n_nopol  = std(nv(t_drop+1:T))
std_rb_nopol = std(rbv(t_drop+1:T))
std_pr_nopol = std(prv(t_drop+1:T))
std_pib_nopol= std(pibv(t_drop+1:T))
std_d_nopol  = std(dv(t_drop+1:T))

%% Pick up bank run events

window = 8;  % how many periods before and after runs to plot

% threshold net worth to define 'normal' times
%thre_n = n_ss;
thre_n = mean(nv(1000:T)) - 2*std(nv(1000:T));

% pick up all runV = 1

num_norun  = 0;
num_run    = 0;
num_nega   = 0;

normal_num = 0;
total_num  = 0;

sum_l_norun  = 0;
sum_l_run    = 0;
sum_l_normal = 0;
sum_p_normal = 0;

rkhat_run  = zeros(1);
rkhats_run = zeros(1);
h_run   = zeros(1);
y_run   = zeros(1);
c_run   = zeros(1);
i_run   = zeros(1);
k_run   = zeros(1);
pib_run = zeros(1);
n_run   = zeros(1);
l_run   = zeros(1);
rb_run  = zeros(1);
pr_run  = zeros(1);
shock_run = zeros(1);
xi_run  = zeros(1);
tfp_run = zeros(1);
pibv_nega = zeros(1);
nv_nega   = zeros(1);

for t = t_drop:T-window

    if xiv(t) < 1
        total_num = total_num + 1;
    end

    if xiv(t)<1 %&& nv(t-1) > thre_n

        num_run = num_run + 1;

        rkhat_run(num_run,1:1+2*window)  = rkhat(t-window:t+window);
        rkhats_run(num_run,1:1+2*window) = rkhats(t-window:t+window);
        h_run(num_run,1:1+2*window)   = hv(t-window:t+window);
        y_run(num_run,1:1+2*window)   = yv(t-window:t+window);
        c_run(num_run,1:1+2*window)   = cv(t-window:t+window);
        i_run(num_run,1:1+2*window)   = iv(t-window:t+window);
        k_run(num_run,1:1+2*window)   = kv(t-window:t+window);
        pib_run(num_run,1:1+2*window) = pibv(t-window:t+window);
        n_run(num_run,1:1+2*window)   = nv(t-window:t+window);
        l_run(num_run,1:1+2*window)   = lv(t-window:t+window);
        rb_run(num_run,1:1+2*window)  = rbv(t-window:t+window);
        pr_run(num_run,1:1+2*window)  = prv(t-window:t+window);
        shock_run(num_run,1:1+2*window) = shocks(t-window:t+window);
        tfp_run(num_run,1:1+2*window) = av(t-window:t+window);
        xi_run(num_run,1:1+2*window)  = xiv(t-window:t+window);

        sum_l_run   = sum_l_run + lv(t);

    else

        sum_l_norun = sum_l_norun + lv(t);

        if pibv(t) > 0

            num_norun   = num_norun + 1;
        else

            num_nega   = num_nega + 1;
            pibv_nega(num_nega) = pibv(t);
            nv_nega(num_nega)   = nv(t);

        end
    end

%    if lv(t) < thre_l
    if nv(t) > thre_n
        normal_num   = normal_num + 1;
        sum_l_normal = sum_l_normal + lv(t);
        sum_p_normal = sum_p_normal + prv(t);
    end

end

uncon_run_all = total_num / (T-window-t_drop)
avg_l_run     = sum_l_run   / num_run;
avg_l_norun   = sum_l_norun / num_norun;
avg_l_all     = mean(lv(1000:T));
avg_l_normal  = sum_l_normal / normal_num;
avg_pr_normal = sum_p_normal / normal_num;

run_avg_rkhat  = mean(rkhat_run,1);
run_avg_rkhats = mean(rkhats_run,1);
run_avg_h   = mean(h_run,1);
run_avg_y   = mean(y_run,1);
run_avg_c   = mean(c_run,1);
run_avg_i   = mean(i_run,1);
run_avg_k   = mean(k_run,1);
run_avg_pib = mean(pib_run,1);
run_avg_n   = mean(n_run,1);
run_avg_l   = mean(l_run,1);
run_avg_rb  = mean(rb_run,1);
run_avg_pr  = mean(pr_run,1);
run_avg_shock = mean(shock_run,1);
run_avg_xi  = mean(xi_run,1);
run_avg_tfp = mean(tfp_run,1);

%% pick up runV = 1 after no run some periods

norun_period = 12;
run_before   = 8;
run_after    = 20;

run_num2 = 0;

rkhat_run2  = zeros(1);
rkhats_run2 = zeros(1);
h_run2   = zeros(1);
y_run2   = zeros(1);
c_run2   = zeros(1);
i_run2   = zeros(1);
k_run2   = zeros(1);
pib_run2 = zeros(1);
n_run2   = zeros(1);
l_run2   = zeros(1);
rb_run2  = zeros(1);
pr_run2  = zeros(1);
shock_run2 = zeros(1);
tfp_run2 = zeros(1);
xi_run2  = zeros(1);
d_run2   = zeros(1);
spread_run2 = zeros(1);
erk_run2    = zeros(1);

norm_num  = 0;
norm_nega = 0;
norm_posi = 0;
norm_y   = zeros(1,1);
norm_c   = zeros(1,1);
norm_l   = zeros(1,1);
norm_k   = zeros(1,1);
norm_i   = zeros(1,1);
norm_n   = zeros(1,1);
norm_pr  = zeros(1,1);
norm_rb  = zeros(1,1);
norm_h   = zeros(1,1);
norm_pib = zeros(1,1);
norm_d   = zeros(1,1);
norm_rkhat   = zeros(1,1);

for t = t_drop:T-run_after

    if sum(runv(t-norun_period:t-1)) == 0
        
        if xiv(t) < 1

        run_num2 = run_num2 + 1;

        rkhat_run2(run_num2,1:1+run_before+run_after)  = rkhat(t-run_before:t+run_after);
        rkhats_run2(run_num2,1:1+run_before+run_after) = rkhats(t-run_before:t+run_after);
        h_run2(run_num2,1:1+run_before+run_after)      = hv(t-run_before:t+run_after);
        y_run2(run_num2,1:1+run_before+run_after)      = yv(t-run_before:t+run_after);
        c_run2(run_num2,1:1+run_before+run_after)      = cv(t-run_before:t+run_after);
        i_run2(run_num2,1:1+run_before+run_after)      = iv(t-run_before:t+run_after);
        k_run2(run_num2,1:1+run_before+run_after)      = kv(t-run_before:t+run_after);
        pib_run2(run_num2,1:1+run_before+run_after)    = pibv(t-run_before:t+run_after);
        n_run2(run_num2,1:1+run_before+run_after)      = nv(t-run_before:t+run_after);
        l_run2(run_num2,1:1+run_before+run_after)      = lv(t-run_before:t+run_after);
        rb_run2(run_num2,1:1+run_before+run_after)     = rbv(t-run_before:t+run_after);
        pr_run2(run_num2,1:1+run_before+run_after)     = prv(t-run_before:t+run_after);
        shock_run2(run_num2,1:1+run_before+run_after)  = shocks(t-run_before:t+run_after);
        tfp_run2(run_num2,1:1+run_before+run_after)    = av(t-run_before:t+run_after);
        xi_run2(run_num2,1:1+run_before+run_after)     = xiv(t-run_before:t+run_after);
        d_run2(run_num2,1:1+run_before+run_after)      = dv(t-run_before:t+run_after);
%        spread_run2(run_num2,1:1+run_before+run_after) = exp(rkhat(t-run_before:t+run_after))-rbv(t-run_before:t+run_after);
        spread_run2(run_num2,1:1+run_before+run_after) = exp(erk(t-run_before:t+run_after))-rbv(t-run_before:t+run_after);
        erk_run2(run_num2,1:1+run_before+run_after)    = erk(t-run_before:t+run_after);
        
        else
            
        norm_num = norm_num + 1;
        norm_y(norm_num)   = yv(t);
        norm_c(norm_num)   = cv(t);
        norm_l(norm_num)   = lv(t);
        norm_k(norm_num)   = kv(t);
        norm_i(norm_num)   = iv(t);
        norm_n(norm_num)   = nv(t);
        norm_pr(norm_num)  = prv(t);
        norm_rb(norm_num)  = rbv(t);
        norm_h(norm_num)   = hv(t);
        norm_pib(norm_num) = pibv(t);
        norm_d(norm_num)   = dv(t);
        norm_rkhat(norm_num)   = rkhat(t);

        if pibv(t) < 0
            norm_nega = norm_nega+1;
        else
            norm_posi = norm_posi + 1;
        end
        
        end
    end
end

run_norun_before = run_num2 / (T-window-t_drop)

run2_avg_rkhat  = mean(rkhat_run2,1);
run2_avg_rkhats = mean(rkhats_run2,1);
run2_avg_h   = mean(h_run2,1);
run2_avg_y   = mean(y_run2,1);
run2_avg_c   = mean(c_run2,1);
run2_avg_i   = mean(i_run2,1);
run2_avg_k   = mean(k_run2,1);
run2_avg_pib = mean(pib_run2,1);
run2_avg_n   = mean(n_run2,1);
run2_avg_l   = mean(l_run2,1);
run2_avg_rb  = mean(rb_run2,1);
run2_avg_pr  = mean(pr_run2,1);
run2_avg_shock = mean(shock_run2,1);
run2_avg_tfp = mean(tfp_run2,1);
run2_avg_xi  = mean(xi_run2,1);
run2_avg_d   = mean(d_run2,1);
run2_avg_spread = mean(spread_run2,1);
run2_avg_erk = mean(erk_run2,1);

norm_y_avg   = mean(norm_y);
norm_c_avg   = mean(norm_c);
norm_l_avg   = mean(norm_l);
norm_k_avg   = mean(norm_k);
norm_i_avg   = mean(norm_i);
norm_n_avg   = mean(norm_n);
norm_pr_avg  = mean(norm_pr);
norm_rb_avg  = mean(norm_rb);
norm_h_avg   = mean(norm_h);
norm_pib_avg = mean(norm_pib);
norm_d_avg   = mean(norm_d);
norm_rkhat_avg   = mean(norm_rkhat);

nega_prob = norm_nega/norm_num;

%% take percentile around the mean path

range = 10;

y_b  = prctile(y_run2,range,[1]);
y_a  = prctile(y_run2,100-range,[1]);
c_b  = prctile(c_run2,range,[1]);
c_a  = prctile(c_run2,100-range,[1]);
i_b  = prctile(i_run2,range,[1]);
i_a  = prctile(i_run2,100-range,[1]);
h_b  = prctile(h_run2,range,[1]);
h_a  = prctile(h_run2,100-range,[1]);
l_b  = prctile(l_run2,range,[1]);
l_a  = prctile(l_run2,100-range,[1]);
n_b  = prctile(n_run2,range,[1]);
n_a  = prctile(n_run2,100-range,[1]);
pr_b = prctile(pr_run2,range,[1]);
pr_a = prctile(pr_run2,100-range,[1]);
rb_b = prctile(rb_run2,range,[1]);
rb_a = prctile(rb_run2,100-range,[1]);
tfp_b = prctile(tfp_run2,range,[1]);
tfp_a = prctile(tfp_run2,100-range,[1]);
d_b  = prctile(d_run2,range,[1]);
d_a  = prctile(d_run2,100-range,[1]);
spread_b  = prctile(spread_run2,range,[1]);
spread_a  = prctile(spread_run2,100-range,[1]);
erk_b  = prctile(erk_run2,range,[1]);
erk_a  = prctile(erk_run2,100-range,[1]);

norm_pib_over_k = mean(norm_pib./norm_k);
norm_pib_over_n = mean(norm_pib./norm_n);

%% plot histogram for pr and pd

close all

figure('Units','Inches','OuterPosition',[1 0 7 7])
hist  = histogram(pr_next(1000:end),100,'Normalization','Probability','FaceAlpha',0.5,'LineWidth',1,'EdgeAlpha',0);
%hist  = histogram(pr_next(1000:end),(0:0.005:0.5),'Normalization','Probability','FaceAlpha',0.5,'LineWidth',1,'EdgeAlpha',0);
%ylim([0 0.9]);
%xlim([0 0.5]);
%xticks(0:0.1:0.5)
grid on
set(gca,'FontSize',30);
set(gca,'yticklabel',num2str(get(gca,'ytick')','%4g'))

figure('Units','Inches','OuterPosition',[1 0 7 7])
hist  = histogram(pd_next(1000:end),100,'Normalization','Probability','FaceAlpha',0.5,'LineWidth',1,'EdgeAlpha',0);
%hist  = histogram(pd_next(1000:end),(0.5:0.005:1),'Normalization','Probability','FaceAlpha',0.5,'LineWidth',1,'EdgeAlpha',0);
%ylim([0 0.9]);
%xlim([0.5 1.0]);
%xticks(0.5:0.1:1)
grid on
set(gca,'FontSize',30);
set(gca,'yticklabel',num2str(get(gca,'ytick')','%4g'))

%% plot histogram for insolvency

figure('Units','Inches','OuterPosition',[1 0 7 7])
hist  = histogram(exp(rkhat(1000:end))-exp(rkhat_insol(1000:end)),100,'Normalization','Probability','FaceAlpha',0.5,'LineWidth',1,'EdgeAlpha',0);
%ylim([0 0.16]);
grid on
set(gca,'FontSize',30);
set(gca,'yticklabel',num2str(get(gca,'ytick')','%4g'))

return

%% clear large data

clear shocks kv rbv nv xiv hv cv yv iv lv pibv runv rhseev rkhat rkhats prv av erk x x_nonlog x_allnonlog dv
clear rkhat_run rkhats_run h_run y_run c_run i_run k_run pib_run n_run l_run rb_run pr_run shock_run xi_run tfp_run pibv_nega nv_nega
clear rkhat_run2 rkhats_run2 h_run2 y_run2 c_run2 i_run2 k_run2 pib_run2 n_run2 l_run2 rb_run2 pr_run2 shock_run2 tfp_run2 xi_run2 d_run2
clear norm_y norm_c norm_l norm_k norm_i norm_n norm_pr norm_rb norm_h norm_pib norm_l_posi norm_d norm_rkhat
